A projection method for statics and dynamics of lattice spin systems 
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A method based on Monte Carlo sampling of the probability flows projected onto the subspace of 
one or more slow variables is proposed for investigation of dynamic and static properties of lattice 
spin systems. We illustrate the method by applying it, with projection onto the order-parameter 
subspace, to the three-dimensional 3-state Potts model in equilibrium and to metastable decay in a 
three-dimensional 3-state kinetic Potts model. 
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The widespread use of computer simulations in many 
fields of physics presents a constant challenge to develop 
new, faster, and more efficient algorithms. Here we de- 
scribe a method to study dynamic and static properties 
of spin lattice systems. It is based on Monte Carlo sam- 
pling of the probability flow projected onto the subspace 
of one or more slow variables. Here we use the order pa- 
rameter. The projected information is subsequently used 
to reconstruct dynamic and/or static quantities. This 
idea was first explored by Schulman and a similar 
approach for the dynamics of metastable decay was de- 
veloped by Lee et al. and in Ref . H . The present work 
extends these developments. We show how appropriately 
sampled projected probability flows can be utilized to in- 
vestigate static probability distributions, as well as the 
dynamics of spin systems. The method is illustrated for 
three cases. The first two deal with three-dimensional, 
three-state Potts models in equilibrium. The ferromag- 
netic model at its transition temperature is considered in 
order to explain, for a generic situation, the basics of the 
method. This model is of interest also in lattice-gauge 
theory Q. Next, the antiferromagnetic model below its 
critical temperature is included to show an unusual ex- 
ample of the phase structure. Although only intended 
as an illustration, these results greatly elucidate previ- 
ous Monte Carlo observations of a medium-temperature 
phase in this model ■ Our third example concerns the 
application to the dynamics of metastable decay in the 
3-dimensional 3-state kinetic ferromagnetic Potts model, 
which can be regarded as an extension of the Ising model 
approximation for extremely anisotropic magnetic sys- 
tems . We demonstrate the strength of the method 
by measuring metastable lifetimes in a region of weak 
fields, where direct simulations are not feasible. 

Consider a 3-state Potts model with the Hamiltonian 
H = —JJ2{i j) &{ a U a j) where Oi £ {0, 1, 2} is the "spin" 
at lattice site i, and the summation runs over all nearest- 
neighbor pairs on a simple-cubic lattice. Macroscopi- 
cally, the system is characterized by the concentrations 
{no, n 2 }, X) n i— 1) °f spins in the three states. These 



triples are mapped into an equilateral triangle represent- 
ing the order-parameter space. In Figs. 1 and 2, the 
projection of a point onto the axis extending from the 
j-th corner gives the concentration rn. 

Consider a simulation using a Monte Carlo method 
with local updates at randomly chosen sites. At any given 
moment, the spins can be divided into classes specified by 
the state a of the spin and by the numbers {a, b, 6 — a — b} 
of its neighbors in the states {0, 1, 2}. Let c° b (no, n\, n 2 ) 
be the average equilibrium population of spins in the class 
{(j, a, b}, conditional on the total concentrations rii of 
spins in state i. Further, denote by the probability 
that a spin in the class {a, a, 6} will flip to the state o' 
when visited by the updating algorithm. 

The central objects of the proposed method are the 
global flip rates v acr >(no, ni , 712): 

*w(no,ni,n 2 ) = ^c£ 6 (no,ni,n 2 ) Kb ■ (!) 

a b 

They define an ideal aggregate Markov process Q on 
the order-parameter space {(n Q ,ni,n 2 )}. This pro- 
cess gives an exact equilibrium aggregated probability 
density P(tiq, n\, n 2 ) ||, which can be calculated from 
the detailed-balance condition P(riQ + 1, n\, ^2)^01(^0 + 
1, ni, ^2) = P{no, ni + 1, n 2 )v w (n , ni + 1, n 2 ) plus anal- 
ogous equations for the other transitions JTo| | . 

Of special importance are the zeroes of the "drift func- 
tions" v aa i — v a i a , since they determine the extremal 
points of the probability densities. For the statics of the 
Potts model, we restrict ourselves to these zero loci. 

We simulated the ferromagnetic model (J=l) at its 
first-order phase transition temperature T=l. 81618 |ll]] 
for a 20 3 lattice. Flip-rate histograms were binned us- 
ing 2 14 equilateral triangles covering the whole order- 
parameter triangle, and 10 6 configurations were typically 
generated. Figure [l](a) shows a schematic picture of the 
drift in one of the three directions. Its zero locus consists 
of two parts. The first one is the symmetry axis of the tri- 
angle, and the second is a symmetric arc. We only show 
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a portion of the arc, since the statistics become insuf- 
ficient far from the probability density maxima. Flows 
in different directions are related by rotations of 2ir/3. 
All three zero loci are superimposed in Fig. |l|(b). Where 
three loci corresponding to different directions intersect, 
the probability density has a stationary point. There 
are three stable extrema Si, corresponding to the three 
ordered phases. The three points Ui are saddle points, 
and the center of the triangle represents the stable, dis- 
ordered phase. This example shows how one identifies 
candidates for stable phases from the intersections of the 
drift-function zero loci. 




FIG. 1. (a): Zeros of the drift V12 — V21 in the 1 <-> 2 direc- 
tion for the 3-dimensional 3-state Potts ferromagnet. Light 
solid lines are stable loci, while dashed lines are unstable. 
Heavy arrows indicate the direction of the probability flow, 
(b): Zero loci for the three directions superimposed. See text 
for a full description. 

Next, we turn to the antiferromagnetic Potts model 
(J=— 1). Because of the sublattice-symmetry breaking, 
we sampled the flip rates for each sublattice separately. 
Figure || shows the union of the zero loci of the drift 
functions in all three directions. As in the ferromagnetic 
case, for each direction there is the symmetry induced 
straight line (triangle axis) and a nontrivial part, which 
is here a closed curve. The closed-curve parts of the zero 
loci have several interesting properties: 

1. The closed loci for the three different directions are 
identical and can be accurately parameterized by 

{m cos t + r cos 2t, m sin t — r sin 2t} , t £ (0, 2ir) . (2) 

2. Positions of the sublattices on this curve are corre- 
lated in such a way that their distance (sublattice mag- 
netization difference \AB\ — 2m in Fig. ||) is constant. 

3. There are finite-size effects in the diameter and 
shape of the curve, but there is no sign that it separates 
into three distinct components. 

We stress that the only observed deviations from these 
properties are numerical uncertainties on the order of 



the discreteness of the order-parameter space. On lat- 
tices smaller than 32 3 , we sampled complete flip-rate his- 
tograms as in the ferromagnetic case. On larger lattices, 
we only measured flip rates at isolated points of the order- 
parameter space (typically 5 x 10 4 configurations). The 
data were then interpolated and used to find the roots 
of the drift functions. In this way, we could confirm the 
above observations, even with lattices as large as 64 3 , by 
measuring flip rates only in the vicinity of the expected 
zero locus. In fact, having located two points we can 
predict the rest with an accuracy suggesting that the pa- 
rameterization of Eq. (0) may be exact. 




FIG. 2. Zero locus of the drift functions v aa i — v a i a of 
the 3-dimensional 3-state Potts antiferromagnet at T = 1.0 J. 
The dashed sections of the triangle axes are unstable in direc- 
tions perpendicular to the axes. The solid parts are stable. 
The closed, solid curve is stable in all three directions. This 
curve, obtained for a 64 3 lattice, represents the degenerate 
maximum of the probability density of a sublattice location 
in the order-parameter space. 

The above observations mean that the probability den- 
sity has a degenerate maximum. Whereas the magnitude 
of the sublattice magnetization difference is fixed, it can 
point in an arbitrary direction. Thus, a system with only 
three microscopic states exhibits a continuous symme- 
try on the macroscopic scale. This is the rotationally 
symmetric phase which was found in recent Monte Carlo 
simulations However, the conventional probability 

density sampling technique did not provide completely 
convincing evidence that this symmetry is not (weakly) 
broken. The present approach gives much stronger sup- 
port for the symmetric phase. It locates the extrema 
of the probability density without sampling the distribu- 
tion itself. With relatively modest statistics, the sym- 
metric property was corroborated with an accuracy ap- 
proaching the limit imposed by the discreteness of the 
order-parameter space. Such precision would require an 
enormous effort if one were to investigate the sampled 
distribution directly. 

Our final example deals with metastable decay in the 
ferromagnetic Potts model with the Hamiltonian H = 
- J T,(i,j) S (^i^j) + H T,i[ s (°^j) -H^Vj)]- Here, we 
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have added a term describing the interaction with the ex- 
ternal field H . With this choice of the external field, the 
model can be regarded as an extension of the kinetic Ising 
model, which can serve as an approximation for nanoscale 
ferromagnets The third spin state of the present 

model can mimic local magnetization "perpendicular" to 
the external field and allows for "finite anisotropy." The 
Glauber dynamics with updates at randomly chosen sites 
is used throughout the rest of the paper. Time is mea- 
sured in Monte Carlo Steps per Spin (MCSS). 

An essential quantity related to metastability is the 
lifetime r, defined as follows. All spins are initialized 
in state 0, the temperature T < T c is fixed, and an ex- 
ternal magnetic field H favoring state f and disfavoring 
state is applied. The field does not interact with spins 
in state 2. This initial state is metastable and decays 
through the nucleation and subsequent growth of stable- 
phase droplets. The average time needed to reach a con- 
figuration with half of the system in the stable phase is r. 
The difficulty is that realistic models, subjected to "ex- 
perimentally reasonable" magnetic fields, have lifetimes 
which are extremely long in terms of the Monte Carlo 
time. Here we describe a significant extension of the 
method proposed in Refs. jj]-[3). This allows us to obtain 
lifetimes in arbitrarily weak fields without prohibitively 
lengthy simulations. 

The projected flip rates are defined, as in the static 
case, in terms of the normalized spin-class populations 
c ab( n ) an d flipping probabilities : 

9in) = £ cl b {n) p% , s(n) = £ c l ab {n) p%. (3) 

ab,a^l ab, cr^l 

We parameterize g and s only by the total number n 
of spins in state f, thus projecting the population data 
onto a one-dimensional histogram. The rates g and s 
correspond to growth and shrinkage of the stable phase. 
They depend on the external field and on the way the 
configurations are generated, as explained below. 

The flip rates are used to map the metastable-decay 
dynamics onto a one-dimensional absorbing Markov 
chain. We assign to all configurations with n overturned 
spins a single state n in the chain. The one-dimensional 
dynamics is given by the flip rates. From state n we 
have the probability g(n) of jumping to state n + 1, the 
probability s(n) of jumping to n — 1, and the probability 
l — s(n)—g(n) of remaining in the current state. This ran- 
dom walk starts at n — and terminates when it reaches 
n = N, corresponding to a stopping criterion which we 
chose to be N = V/2, with V the volume of the system. 
Using standard methods from the theory of absorbing 
Markov chains |^,|f^| , we obtain the mean lifetime r and 
the total average time h(n) (measured in MCSS, with 
h(N) = 0) spent by the random walker in the state n as 

E 1 , f \ ii i\ V^ 1 + s(n)h(n) 
h{n) , h(n-l) = y _> . (4) 

n=0 y ^ ' 



How accurately these formulas reproduce the lifetime de- 
pends on how the class populations are sampled. One op- 
tion is to sample them in zero external field in an equilib- 
rium ensemble with conserved order parameter for each 
needed value of n plpl. Such data can be used to esti- 
mate the lifetime in very weak fields, but in strong fields 
it is underestimated because the class populations near 
the top of the free-energy barrier are not reproduced well 
by the equilibrium ensemble. 

The simple but important improvement presented here 
is the way the class populations are measured. At any 
time, the system is only allowed to have n, the num- 
ber of spins in the stable phase, larger than a time- 
dependent lower bound, n m j n . Simulation starts with 
n m in(t=0) = — 1, and n m ; n is increased slowly. The class 
populations c G ah are sampled during this forced-escape 
simulation with an applied external field and are sub- 
sequently used to calculate the lifetime from Eqs. (||JJ). 
Why does this work? Consider first the limit of zero 
forcing speed (n m m=— f at all times). One starts a sim- 
ulation by releasing a "random walker" from the initial 
state. When the walker reaches the absorbing state, one 
starts a new run and repeats the whole process until 
iVe S c escapes from metastability are realized. One can 
imagine that each walker generates an oriented "world 
line." Therefore, an "equation of continuity" holds for 
each n < N: N n ^ n+1 = N csc + N n+1 ^ n with Ni->j the 
number of transitions between the subspaces with i and 
j overturned spins, ft is straightforward to express this 
equality in terms of the transition probabilities and class 
populations as generated by the simulation. One obtains 
a relation for h(n) equivalent to Eq. (^) |l(J. Thus, if the 
forcing is infinitely slow, Eqs. (0,0) give the exact r. 

Now, if the rate of increase of n m i n is nonzero but 
sufficiently small, the system always produces nearly 
the "correct" configurations as if there were no forcing. 
While deep in the metastable free-energy well, forcing 
prevents the system from returning to n < n m i n , but it 
still allows it to thermalize and generate metastable con- 
figurations with n at and slightly above n min . As n min in- 
creases, the procedure scans the configurations along the 
escape path from metastability. When n min approaches 
the top of the free-energy barrier, the system has a better 
chance to escape, and the stable phase grows too quickly 
to equilibrate. There is nothing to prevent escape, be- 
cause there is no upper bound on n which would cause 
unwanted thermalization. The system is free to escape 
through natural noncquilibrium configurations. 

To approach the slow-forcing limit in practice, one per- 
forms a series of measurements and determines a value of 
dn m i n /dt below which the estimated lifetime is insensitive 
to the choice of dn m i n /dt. We expect such forcing speed 
to be related to the inverse equilibration time for states 
deep in the metastable well. Another computationally 
important aspect is to obtain sufficient statistics for the 
forced escapes, since the class-population data tend to 
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be noisy at and beyond the top of the free-energy bar- 
rier. Here, we measured about 10 3 escapes at rates up to 




FIG. 3. The metastable lifetime of a 3-dimensional 3-state 
kinetic Potts ferromagnet as a function of the magnetic field 
for two system sizes at T ~ 0.55T C . Symbols o and □ are 
simulation results (with error bars smaller than the symbol 
size), and lines connect the points (+) calculated from the 
class-population data sampled by the forced-escape method. 
The inset shows a comparison of the measured (circles) and 
calculated (line) relative standard deviation of the lifetime. 

Figure || shows a comparison between the lifetimes ob- 
tained from direct simulations and those calculated with 
the forced-escape method. The agreement is excellent for 
lifetimes in the whole region accessible to direct simula- 
tions. Moreover, the forced escape method can provide 
lifetime estimates deep in the region of weak fields, where 
direct simulation is practically impossible. We empha- 
size that having measured the class populations, one can 
utilize them to calculate lifetimes for an arbitrary local 
dynamic with random updates. We can e.g. calculate 
what the lifetimes would be if we used the Metropolis in- 
stead of the Glauber dynamic simply by replacing the flip 
probabilities . Although different dynamics produce 
different flip rates, the class-populations remain close to 
local equilibrium and are therefore similar for all dynam- 
ics that obey detailed balance. 

In a similar way as the average lifetime r, one can use 
the flip rates to calculate higher moments of the lifetime 
probability distribution [0. The inset in Fig. || shows 
the relative standard deviation Ar/r as a function of 
the field. Despite the fact that the higher moments, un- 
like the mean lifetime, are only approximate even in the 
slow- forcing limit, the agreement between the directly 
simulated results and those based on the forced-escape 
method (solid line) is good. We believe this is because the 
relevant Markov chains are nearly weakly lumpable p|Jl3| 
with respect the states along the escape from metastabil- 
ity. For a quantitative description of the form of r and 
Ar/r, see Refs. § and g. 



In summary, we propose a new method to study the 
dynamic and static properties of lattice spin systems. It 
consists in Monte Carlo sampling of the spin-class popu- 
lations in a projected subspace. These are subsequently 
used to calculate the spin- flip rates. In the static case, the 
knowledge of the flip rates is equivalent to the informa- 
tion contained in the probability density. However, the 
accuracy of the method is much better than with direct 
distribution sampling. We have demonstrated this on the 
example of the 3-state antiferromagnetic Potts model. 

The flip rates can be also utilized to obtain dynamic 
characteristics of the systems, such as lifetimes, in very 
weak fields where ordinary simulations are not feasible. 
Although the resulting dynamic is only an approxima- 
tion, it provides accurate estimates. The detailed infor- 
mation about the class populations enables one to calcu- 
late the lifetime for an arbitrary dynamic with updates 
at randomly chosen sites, independently of the dynamic 
used during the data sampling. 
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